% this is a simpler 2-D diffusion problem so that I can
% ensure I get bc's right

% looks like i got the bc's set, but i can easily simplify the code
% things to do:
%   * clean up repetition
%   * variable material stuff, probably copy & paste from 1st run
%   * coarse mesh material assigments with fine mesh intervals defined
%   * multigroup formulation

clear; clc;

% something like water
dc = 0.20;
sa = 0.02;

xlen = 5;  % cm
ylen = 5;  % cm
N = 10; % number of x intervals
M = 10; % number of y intervals
dx = xlen/N;
dy = ylen/M;

% boundary conditions: left, right, bottom, top; 0 = vacuum, 1 = reflect
BCL = 0;
BCR = 0;
BCB = 1;
BCT = 1;


src = zeros(N,M); % uniform volume source
src(1:25,1:25)=1; 

% there are (N+1)*(M+1) fluxes to define
%A = zeros( (N+1)*(M+1) );
lenk = (N+1)*(M+1);
AC = zeros( lenk, 1 );
AL = zeros( lenk-1, 1 );
AR = AL;
AB = zeros( lenk-N-1, 1 );
AT = AB;
S = zeros( (N+1)*(M+1), 1 );

% interior coefficients
for i = 2:N 
    for j = 2:M
        k = i+(j-1)*( N+1 ); % 10 being the number of x ints
        alft = -dc/dx^2;
        argt = -dc/dx^2;
        abot = -dc/dy^2;
        atop = -dc/dy^2;    
        a    = sa - (alft+argt+abot+atop);
        AL( k-1 ) = alft;  
        AR( k ) = argt;   
        AB( k-N-1 ) = abot;    
        AT( k ) = atop;
        AC( k ) = a;
        S(k)    = src(i,j);              
    end
end

% left and right edges
for j = 2:M
    i = 1;
    k = i+(j-1)*( N+1 );
    if BCL == 1
        AR( k ) = -1;             
        AC( k ) = 1;  
    else
        AR( k ) = -0.5*dc/dx;    
        AC( k ) = 0.25+0.5*dc/dx;  
    end
    i = N+1;
    k = i+(j-1)*( N+1 );
    if BCR == 1
        AL( k-1 ) = -1;              
        AC( k ) = 1;
    else
        AL( k-1 ) = -0.5*dc/dx;
        AC( k ) = 0.25+0.5*dc/dx;      
    end
end
% top and bottom edges
for i = 2:N
    j = 1;
    k = i+(j-1)*( N+1 );  
    if BCB == 1
        AT( k )     = -1;             
        AC( k  )    = 1;  
    else
        AT( k )     = -0.5*dc/dy;     
        AC( k )     = 0.25+0.5*dc/dy;        
    end 
    j = M+1;
    k = i+(j-1)*( N+1 );
    if BCT == 1
        AB( k-N-1 ) = -1;            
        AC( k )     = 1; 
    else
        AB( k-N-1 ) = -0.5*dc/dy;   
        AC( k )     = 0.25+0.5*dc/dy;        
    end
end

% finally, the corners
i=1; j=1; k = i+(j-1)*( N+1 );      % BOTTOM LEFT
if (BCL==1&&BCB==1)
    AR( k ) = -1; %right  
    AT( k ) = -1; %top
    AC( k ) = 2;  %center
elseif(BCL==1&&BCB==0)
    AR( k ) = -1; %right  
    AC( k ) = 1; %center   
elseif (BCL==0&&BCB==1)
    AT( k ) = -1; %top
    AC( k ) = 1; %center   
else    
    AR( k ) = -0.5*dc/dx; %right  
    AT( k ) = -0.5*dc/dy; %top
    AC( k ) = 0.5 + 0.5*(dc/dx+dc/dy); %center  
end
i=1; j=M+1; k = i+(j-1)*( N+1 );    % TOP LEFT
if (BCL==1&&BCT==1)
    AR( k )     = -1; %right  
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 2;  %center
elseif(BCL==1&&BCT==0)
    AR( k ) = -1; %right  
    AC( k ) = 1;  %center   
elseif (BCL==0&&BCT==1)
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 1;  %center   
else    
    AR( k )     = -0.5*dc/dx; %right  
    AB( k-N-1 ) = -0.5*dc/dy; %bot
    AC( k )     = 0.5 + 0.5*(dc/dx+dc/dy); %center  
end
i=N+1; j=1;  k = i+(j-1)*( N+1 );   % BOTTOM RIGHT
if (BCR==1&&BCB==1)
    AL( k-1 ) = -1; %left  
    AT( k )    = -1; %top
    AC( k )   = 2;  %center
elseif(BCR==1&&BCB==0)
    AL( k-1 ) = -1; %left  
    AC( k )   = 1; %center   
elseif (BCR==0&&BCB==1)
    AT( k )     = -1; %top
    AC( k )                   = 1; %center   
else    
    AL( k-1 ) = -0.5*dc/dx; %left  
    AT( k )   = -0.5*dc/dy; %top
    AC( k )   = 0.5 + 0.5*(dc/dx+dc/dy); %center  
end
i=N+1; j=M+1;  k = i+(j-1)*( N+1 );   % TOP RIGHT
if (BCR==1&&BCT==1)
    AL( k-1 )   = -1; %right  
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 2;  %center
elseif(BCL==1&&BCT==0)
    AL( k-1 ) = -1; %right  
    AC( k )   = 1; %center   
elseif (BCL==0&&BCT==1)
    AB( k-N-1 ) = -1; %bot
    AC( k )     = 1; %center   
else    
    AL( k-1 )   = -0.5*dc/dx; %right  
    AB( k-N-1 ) = -0.5*dc/dy; %bot
    AC( k )     = 0.5 + 0.5*(dc/dx+dc/dy); %center  
end    

%---establish sparse matrix for solution
% pad end of AB,AL
AB = [AB' zeros(1,N+1)]';
AL = [AL' 0]';
% pad top of AR,AT
AT = [zeros(1,N+1) AT']';
AR = [0 AR']';
% form sparse matrix for solution
A = spdiags([AB AL AC AR AT],[-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
% now solve -- will want to implement an iterative solver, i think
phi = A\S;
% generate indices
for i = 1:N+1
    for j = 1:M+1
        k = i+(j-1)*( N+1 );
        indx(k,1)=i;
        indx(k,2)=j;
    end
end
% generate x and y vals
xint = zeros(N+1,1); 
yint=zeros(M+1,1);
for i = 2:N+1
    xint(i)=xint(i-1)+dx;
end
for i = 2:M+1
    yint(i)=yint(i-1)+dy;
end
phimap = zeros(N+1,M+1);
for k = 1:(N+1)*(M+1)
    i = indx(k,1);
    j = indx(k,2); 
    phimap(i,j)=phi(k);
end

surf(xint,yint,phimap')
xlabel('x')
ylabel('y')

